package Models;
 
public class FFT {
 
   int n, m;
   
   // Lookup tables.  Only need to recompute when size of FFT changes.
   double[] cos;
   double[] sin;
 
   double[] window;
   
   public FFT(int n) {
     this.n = n;
     this.m = (int)(Math.log(n) / Math.log(2));
 
     // Make sure n is a power of 2
     if(n != (1<<m))
       throw new RuntimeException("FFT length must be power of 2");
 
     // precompute tables
     cos = new double[n/2];
     sin = new double[n/2];
 
     for(int i=0; i<n/2; i++) {
       cos[i] = Math.cos(-2*Math.PI*i/n);
       sin[i] = Math.sin(-2*Math.PI*i/n);
     }
 
     makeWindow();
   }
 
   protected void makeWindow() {
     // blackman window:
     // w(n)=0.42-0.5cos{(2*PI*n)/(N-1)}+0.08cos{(4*PI*n)/(N-1)};
     window = new double[n];
     for(int i = 0; i < window.length; i++)
       window[i] = 0.42 - 0.5 * Math.cos(2*Math.PI*i/(n-1)) 
         + 0.08 * Math.cos(4*Math.PI*i/(n-1));
   }
   
   public double getWindowAcf(){
     return( 2.7932 ); // Blackman window only  
   }
   
   protected void makeNoWindow() {
     window = new double[n];
     for(int i = 0; i < window.length; i++)
       window[i] = 1.0;
   }
   
   public double[] getWindow() {
     return window;
   }
 
   /***************************************************************
   * fft.c
   * Douglas L. Jones 
   * University of Illinois at Urbana-Champaign 
   * January 19, 1992 
   * http://cnx.rice.edu/content/m12016/latest/
   * 
   *   fft: in-place radix-2 DIT DFT of a complex input 
   * 
   *   input: 
   * n: length of FFT: must be a power of two 
   * m: n = 2**m 
   *   input/output 
   * x: double array of length n with real part of data 
   * y: double array of length n with imag part of data 
   * 
   *   Permission to copy and use this program is granted 
   *   as long as this header is included. 
   ****************************************************************/
   public void fft(double[] x, double[] y)
   {
     int i,j,k,n1,n2,a;
     double c,s,e,t1,t2;
   
   
     // Bit-reverse
     j = 0;
     n2 = n/2;
     for (i=1; i < n - 1; i++) {
       n1 = n2;
       while ( j >= n1 ) {
         j = j - n1;
         n1 = n1/2;
       }
       j = j + n1;
     
       if (i < j) {
         t1 = x[i];
         x[i] = x[j];
         x[j] = t1;
         t1 = y[i];
         y[i] = y[j];
         y[j] = t1;
       }
     }
 
     // FFT
     n1 = 0;
     n2 = 1;
   
     for (i=0; i < m; i++) {
       n1 = n2;
       n2 = n2 + n2;
       a = 0;
     
       for (j=0; j < n1; j++) {
         c = cos[a];
         s = sin[a];
         a +=  1 << (m-i-1);
 
         for (k=j; k < n; k=k+n2) {
           t1 = c*x[k+n1] - s*y[k+n1];
           t2 = s*x[k+n1] + c*y[k+n1];
           x[k+n1] = x[k] - t1;
           y[k+n1] = y[k] - t2;
           x[k] = x[k] + t1;
           y[k] = y[k] + t2;
         }
       }
     }
   }                          
 

   // Test the FFT to make sure it's working
   // public static void main(String[] args) {
   public static void main_test(){
     int N = 8;
 
     FFT fft = new FFT(N);
 
     double[] window = fft.getWindow();
     double[] re = new double[N];
     double[] im = new double[N];

     // Impulse
     re[0] = 1; im[0] = 0;
     for(int i=1; i<N; i++)
       re[i] = im[i] = 0;
     beforeAfter(fft, re, im);
 
     // Nyquist
     for(int i=0; i<N; i++) {
       re[i] = Math.pow(-1, i);
       im[i] = 0;
     }
     beforeAfter(fft, re, im);
 
     // Single sin
     for(int i=0; i<N; i++) {
       re[i] = Math.cos(2*Math.PI*i / N);
       im[i] = 0;
     }
     beforeAfter(fft, re, im);
 
     // Ramp
     for(int i=0; i<N; i++) {
       re[i] = i;
       im[i] = 0;
     }
     beforeAfter(fft, re, im);
 
     long time = System.currentTimeMillis();
     double iter = 30000;
     for(int i=0; i<iter; i++)
       fft.fft(re,im);
     time = System.currentTimeMillis() - time;
     System.out.println("Averaged " + (time/iter) + "ms per iteration");
   }
 
   protected static void beforeAfter(FFT fft, double[] re, double[] im) {
     System.out.println("Before: ");
     printReIm(re, im);
     fft.fft(re, im);
     System.out.println("After: ");
     printReIm(re, im);
   }
 
   protected static void printReIm(double[] re, double[] im) {
     System.out.print("Re: [");
     for(int i=0; i<re.length; i++)
       System.out.print(((int)(re[i]*1000)/1000.0) + " ");
 
     System.out.print("]\nIm: [");
     for(int i=0; i<im.length; i++)
       System.out.print(((int)(im[i]*1000)/1000.0) + " ");
 
     System.out.println("]");
   }
 }
